# define point-quarter distance sampling method
# 
# Author: guochun
###############################################################################



oneSampling.pQuarter=function(distanceOBJ,populationData){
	# select a transect line if the distanceOBJ do not contain this inf.
	if(length(distanceOBJ@slope)==0){
		distanceOBJ=transectLine(distanceOBJ,populationData)
	}
	#select a random point in this transect line
	#add edge adjust method
	#print(distanceOBJ@xrange)
	
	selectx=sample(seq(distanceOBJ@xrange[1]+distanceOBJ@edge,distanceOBJ@xrange[2]-distanceOBJ@edge,1),1)
	selecty=distanceOBJ@slope*selectx+distanceOBJ@b
	win=owin(c(0,populationData@plotdim[1]),c(0,populationData@plotdim[2]))
	select=ppp(x=selectx,y=selecty,window=win,check=F)
	#browser()
	#divide the population into four quadrat center at the random point
	c1=populationData@x > selectx & populationData@y > selecty
	c2=populationData@x > selectx & populationData@y < selecty
	c3=populationData@x < selectx & populationData@y > selecty
	c4=populationData@x < selectx & populationData@y < selecty
	c1.ppp=ppp(x=populationData@x[c1],y=populationData@y[c1],window=win,marks=which(c1),check=F)
	nd1=nncross(select,c1.ppp)
	if(!is.na(nd1[1,2]))
	    availible=c(c1.ppp[nd1[1,2]]$marks)
	c2.ppp=ppp(x=populationData@x[c2],y=populationData@y[c2],window=win,marks=which(c2),check=F)
	nd2=nncross(select,c2.ppp)
	if(!is.na(nd2[1,2]))
	    availible=c(availible, c2.ppp[nd2[1,2]]$marks)
	
	c3.ppp=ppp(x=populationData@x[c3],y=populationData@y[c3],window=win,marks=which(c3),check=F)
	nd3=nncross(select,c3.ppp)
	if(!is.na(nd3[1,2]))
	    availible=c(availible,c3.ppp[nd3[1,2]]$marks)
	c4.ppp=ppp(x=populationData@x[c4],y=populationData@y[c4],window=win,marks=which(c4),check=F)
	nd4=nncross(select,c4.ppp)
	if(!is.na(nd4[1,2]))
	    availible=c(availible,c4.ppp[nd4[1,2]]$marks)
	
	attr(availible,"d")=data.frame(xfrom=selectx,yfrom=selecty,
			x1=noPoint(c1.ppp,nd1,"x"),y1=noPoint(c1.ppp,nd1,"y"),d1=nd1[1,1],
			x2=noPoint(c2.ppp,nd2,"x"),y2=noPoint(c2.ppp,nd2,"y"),d2=nd2[1,1],
			x3=noPoint(c3.ppp,nd3,"x"),y3=noPoint(c3.ppp,nd3,"y"),d3=nd3[1,1],
			x4=noPoint(c4.ppp,nd4,"x"),y4=noPoint(c4.ppp,nd4,"y"),d4=nd4[1,1])
	attr(availible,"distanceOBJ")=distanceOBJ;
	return(availible)
}

noPoint=function(cppp,nd,xy){
	if(is.na(nd[1,2])){
		return(NA)
	}else if(xy == "x"){
		
		return (cppp[nd[1,2]]$x)
	}else{
		return (cppp[nd[1,2]]$y)
	}
		
}


setSampleResult.pQuarter=function(sampleOBJ,selection,populationData){
	return(appendResult(sampleOBJ,selection,populationData))
}



# quick version of point quarter method
samplingQuick.pQuarter=function(sampleMethod,populationData){
	#browser()
	#多取5倍的空白点，最后在删除重复点后，若有多余，则可以随机选择需要的点数
	spt=round(sampleMethod@time*5)
	#if(spt>=length(populationData@x))
	#	spt=length(populationData@x)
	
	repl=TRUE
	rep.time=1
	while(repl & rep.time<=sampleMethod@stopLimit){
		#need a transect line first
		if(sampleMethod@transectLine){
			dline=transectLine(sampleMethod,populationData)
			#while((dline@d/3 -1) <= spt)
			#	dline=transectLine(sampleMethod,populationData)
			
			xrange=dline@xrange
			yrange=dline@slope*xrange+dline@b
			beta=-atan(dline@slope)
			
			#随机选取一个变化方向
			#beta=runif(1,-3*pi/4,3*pi/4)
			A=polarTranf(populationData@x,populationData@y,beta)
			#随机选取一个直线上的点
			#lpoint=polarTranf(runif(1,sampleMethod@edge,populationData@plotdim[1]-sampleMethod@edge),
			#		runif(1,sampleMethod@edge,populationData@plotdim[2]-sampleMethod@edge),beta)
			
			p1=polarTranf(xrange[1],yrange[1],beta)
			p2=polarTranf(xrange[2],yrange[2],beta)
			#为保证每个点间有一定的空隙，用sample命令,但这样产生的点就不是完全随机
			#x=sample(seq(p1[1],p2[1],2),spt)
			x=runif(spt,p1[1],p2[1])
			y=rep(p1[2],times=spt)
			if(all(min(x) >= A[,1]) | all(max(x) <= A[,1]) | all(min(y) >= A[,2]) | all(max(y) <= A[,2])){
				rep.time=rep.time+1
				repl=TRUE
				next()
			}
		# no transect line needed
		}else{
			x=runif(spt,sampleMethod@edge,populationData@plotdim[1]-sampleMethod@edge)
			y=runif(spt,sampleMethod@edge,populationData@plotdim[2]-sampleMethod@edge)
			A=cbind(populationData@x,populationData@y)
		}
		
		B=matrix(c(x,y),nrow=spt)
		knmax=(sampleMethod@k-1)*20+20
		ki=knnx.index(A,B,k=knmax)
		
		result=matrix(NA,nrow=spt,ncol=4)
		#用于管理需要第几个最近邻体
		needadd=matrix(0,nrow=spt,ncol=4)
		for(i in 1:knmax){
			#当无法找到合适的最近临体时，knnx.index会返回-1
		    #这时只要随便付个值，最终不要那个计算的值就行
			den=which(ki[,i]==-1)
			if(length(den)>0){
				ki[den,i]=1
			}
				
			C=A[ki[,i],]
			img=loca(B[,1],B[,2],C[,1],C[,2])
			addloc=1:spt+spt*(img-1)
			needadd[addloc]=needadd[addloc]+1
			#当result中还空，又有合适的临体出现时，才需要付值
			add = (needadd[addloc]>=sampleMethod@k & is.na(result[addloc]))
			if(length(den)>0)
			   add[den]=FALSE
		   
			if(any(add)){
				#kloc=addloc[needadd]%%spt
				#kloc[kloc==0]=spt
				result[addloc[add]]=ki[add,i]
			}
		}
		#没找到某象限内点的
		nNA=which(apply(result,1,function(x) any(is.na(x))))
		#重复的
        temp=table(result[-nNA,])
        #as.numeric(names(temp)[which(temp>1)])
		nl=match(as.numeric(names(temp)[which(temp>1)]),result)%%spt
		del=unique(c(nNA,nl))
		if((spt-length(del))<sampleMethod@time){
			#如果剩下的还不足够符合条件，则在重做一次
			repl=TRUE
			rep.time=rep.time+1
		}else{
			if(length(del)>0){
				result=result[-del,]
				B=B[-del,]	
			}
			select=sample(1:dim(result)[1],sampleMethod@time)
			result=result[select,]
			B=B[select,]
			
			repl=FALSE
		}
	}
	if(rep.time <= sampleMethod@stopLimit){
		d=numeric()
		for(ii in 1:4){
			temp2=A[result[,ii],]
			d=c(d,sqrt((B[,1]-temp2[,1])^2+(B[,2]-temp2[,2])^2))
		}
		return(d)
	}else{
		return(NA)
	}
	
}

polarTranf=function(x,y,beta){
	r=sqrt(x^2+y^2)
	alpha=asin(y/r)
	x=r*cos(alpha+beta)
	y=r*sin(alpha+beta)
	return(matrix(c(x,y),nrow=length(x)))
}


#fist: 3; second:1; third:0; forth:2

loca=function(x0,y0,x1,y1){
	re=rep(0,length(x0))
	re[x1>x0]=re[x1>x0]+2
	re[y1>y0]=re[y1>y0]+1
	re=re+1
	return(re)
}
